February 6, 2008 22:59 WSPC/INSTRUCTION FILE 068'poeschel 



International Journal of Modern Physics C 
© World Scientific Publishing Company 



GRANULAR GAS COOLING AND RELAXATION TO THE 
STEADY STATE IN REGARD TO THE OVERPOPULATED TAIL 
OF THE VELOCITY DISTRIBUTION 



THORSTEN POSCHEL 

Charite, Augustenburger Platz, 10439 Berlin, Germany 
thorsten.poeschel@charite.de 

NIKOLAI V. BRILLIANTOV 

Institute of Physics, University of Potsdam, Am Neuen Palais 10, 14469 Potsdam, Germany 

nbrillia@agnld. uni-potsdam. de 

ARNO FORMELLA 

Universidad de Vigo, Department of Computer Science, Edificio Politecnico, 
32004 Ourense, Spain 
formella@ei.uvigo.es 

Received February 6, 2008 
Revised Day Month Year 

We present a universal description of the velocity distribution function of granular gases, 
/(d), valid for both, small and intermediate velocities where v is close to the thermal 
velocity and also for large v where the distribution function reveals an exponentially 
decaying tail. By means of large-scale Monte Carlo simulations and by kinetic theory 
we show that the deviation from the Maxwell distribution in the high-energy tail leads 
to small but detectable variation of the cooling coefficient and to extraordinary large 
relaxation time. 
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1. Introduction 

Granular gases are characterized by a certain loss of kinetic energy according to the 
dissipative properties of particle collisions. The particle velocities before a collision, 
t/i/2, and after, v[, 2 , are related by the collision rule 

v l = v i 7T~ ( u 12 ' e ) e 

lie (1) 
v 2 = v 2 H — {vi2 • e ) e 

with v\2 = V\ — V2 and the unit vector e = (fi — r%) / \r\ — | at the moment of 
the collision. 
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In the absence of external excitation, a homogeneously initialized granular gas 
stays homogeneous during the first part of its evolution, called homogeneous cooling 
state (HCS). In later stages of its evolution, hydrodynamics instabilities develop 
leading to vortex formation, that is, spacial correlations of the vectorial velocity 
fieldP and finally cluster formation due to a pressure instability.^ In this paper we 
restrict ourselves to the homogeneous cooling state. It should be mentioned that 
heated granular gases with a Gaussian thermostat are equivalent to the HC£p^, 
therefore, the presented results apply also for such systems. 

According to the steady loss of energy, the velocity distribution function, / (v, r), 
depends on time r (measured in units of collisions per particle). Its shape can be 
described by the time-independent function /(c) via^ 

= ^/(*) ( 2 ) 

where 



and v t {t) = v / 2T(r) . (3) 



vt(t) 

The evolution of the granular temperature T(r) is given by Haff 's law^ (with particle 
mass mi = 1), 

J/Tl 

— = -2 7 T, i.e., T(r)=T(0)exp(-2 7 r). (4) 
Here 7 is the cooling coefficient, addressed in detail in the next section. 

2. Universal velocity distribution function 

For small and intermediate velocities, |u| ~ vt, that is, \c\ ~ 1, the scaled velocity 
distribution function /(c) is close to a Gaussian with small deviations that can 
be characterized by the first non-trivial coeffici ent. 02, of a Sonine polynomials 
expansion around the leading Gauss distribution^^ 

/(c) = 7r- 3 / 2 exp(- C 2 )[l + a2 ^ 2 (c 2 )] with S 2 (c 2 ) = \^ - ^ + H (5) 
and 

_ I6(f- £ )(l-2 £ 2 ) 
a2(£) -8I-I7e + 30 £ + 30 2 (l- £ )- (6) 
For large kinetic energy, c> 1, the distribution function decays exponentially 
slow, that is, the distribution function is overpopulated with respect to the Gauss 
distribution one would expect for a molecular gasj^El 

37T 

/(c) = Be~ bc with b = — (7) 

with the second moment of the dimensionless collision integral 

3 

I H a 2 

16 



= - J dc lC \l (/, /) = yfa (1 - £ 2 ) 



(8) 
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and with 

/(/./) = / <lr'2 



de@(-ci2 ■ e ) \ci2 ■ e\ 



1 i 



-3/ (#) /(#)-/(*)/ (30 



(9) 



where c", 2 are the pre-collisional scaled velocities. 

There is a transition region between the ranges of validity of the Sonine expan- 
sion (c « 1) and the expression for the tail (c> 1) where the distribution function 
was quantified in terms of complicated implicit expressions.^ For practical compu- 
tations, however, a much simpler approach may be sufficient: It was shown recently 
that the simple combination of the shown Sonine expansion and the exponential 
tail, 

' Ac 2 e~ c2 [l + a 2 S 2 (c 2 )] for c < c* 
Bc 2 e- bc for c > c* 



f(c) 



(10) 



with the parameters A, B, and the transition velocity c* agrees with large-scale 
DSMC simulation up to a very good accuracy in the entire range of velocities.^ The 
parameters A, B, and c* follow from the conditions of smoothness 

/(c + o) = /(c, - o) 



c,+0 



4/ 
dc 



c,-0 



and normalization: 

b a 2 (2cl - 5c*) 
2 + 2 [1 + a 2 S 2 (cj) 

1 _ fc(c*) , 
1 = 



/J 3 

B = Afc(c*) 



[2 + 6c* (2 + 6c*)]e 



-be 



with 



fc(c*) 



.^Erf(c,)-|[4+a 2 c^(2c^-5)]e- 



[1 + a 2 S 2 (cj)] 



(11) 

(12) 

(13) 
(14) 

(15) 



Solving the fifth order equation (|12p numerically for c*, we obtain A, fc and finally 
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3. Impact of the high-velocity tail on the coefficient of cooling 

By means of the distribution function, Eq. (fTQ| . and its parameters A, B, c* depend- 
ing on the coefficient of restitution e, we can quantify the impact of the exponential 
tail on the cooling coefficient 7, which is one of the most important characteristics 
of a granular gas. Using the standard analysis one can calculate the cooling coef- 
ficient for the homogeneous cooling state, when the stationary velocity distribution 
is achieved: 

J 2 



7=(l-e 2 



12J ' 



(16) 
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where the coefficients Jk depend on the velocity distribution function, 

,h = Jdc lC lc 2 J de6(-c 12 • e)|c 12 • e | fc+1 /(c 1 )/(c 2 ) . (17) 

If we neglect the exponential tail, we obtain the energy decay rate^ 

„(,)_ 1^5* ! + ( 18 ) 

6 i__ a2 ( £ ) 

Complementary, we determin e the c ooling coefficient numerically, using Direct 
Simulation Monte Carlo We initialized a granular gas of N = 10 

particles with velocities drawn from a Gauss distribution at T(0) = 1 and simulated 
until the particle velocities reached the limit of the double precision number repre- 
sentation, i.e., until T « 10~ 23 . For e — 0.9 this corresponds to a total of 5 x 10 10 
collisions or 1,000 collisions per particle. 

For different e we recorded the temperature Tdsmc( t )- Then 7dsmc was 
determined by fitting Tdsmc(t) for r > 1 to its asymptotic law, Tdsmc oc 
exp(-2 7DSMCr), Eq. ©. 

The dependence of 70 on the coefficient of restitution e is shown in Fig.[T]together 
with the numerical result. The numerical results (see below) and the theoretical 
results disregarding the tail, Eq. (fT8|) . virtually coincide. When we plot the difference 
of both curves, however, we find a systematic deviation between the numerical values 
and the analytical result (Fig. [2|) . 

To estimate this rather small difference quantitatively, we write the distribution 
function as a sum of two terms, 

/(c) = /o(c) + /t(c), (19) 




_i 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 



Fig. 1. The cooling coefficient 7(e) as a function of the coefficient of restitution as obtained by 
Direct Simulation Monte Carlo (DSMC) (points) and due to Eq. U8I I. 70. 
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Fig. 2. Difference of the cooling coefficient due to Eq. {18) and the numerically values, (70 — 
7dsmc)i over the coefficient of restitution. Since Eq. (1 1 8 I t disregards the exponential tail of the 
distribution, the shown difference is due to the presence of the tail. The dashed line shows the 
best exponential fit. 



that is, the main part /o(c) (with c extended to infinity), and the pure overpopula- 
tion of the tail, ft(c), 

f (c) = Ac 2 e- c2 [l + a 2 S 2 (c 2 )] (20) 

ft(c)= [Bc 2 exp(-& C )-/o(c)l e(c-c„) . (21) 



The product /(ci)/(c2) in Eq. (Q~7]) reads then 

/(ci)/(c 2 ) = /o(ct)/o(c 2 ) + /o( Cl )/ t (c 2 ) + /t(ci)/o(ca) + /t( Cl )/ t (c 2 ) . 



(22) 



For large c* one can neglect the last term in the above equation. Moreover, in this 
case we approximate ft{c) « £>exp(— be) and c\i ■ e « c\ ■ e, assuming c\ S> c 2 . A 
similar approximation may be applied for the opposite case, c 2 3> c\. Taking then 
into account the symmetry of the integrand in Eq. \Y1\ with respect to c\ and c 2 , 
we finally obtain an estimate for J^: 



71 -2 r(0) 

16" k 



J k = —A* J, 
-A 2 J, 



dc 2 / (c 2 ) 



rfci / de 0(-ci • e ) |ci • e\ k+1 / t (ci 



16 



where is the value of the coefficient Jk for the distribution function /o(c) with 
A in Eq. ([2Qj) equal to 4/y / 7r, that is, for the case when the exponential tail is 
disregarded. In particular, we need 

1 



2^47T 

4 



Be~ bci dcxl-K I sin6Hcos 



ifc+i 



(23) 



J, 



(0) 



7(0) 



2V2tt (^1 
4v / 2^l 1 



16 
3 
1G 



a 2 



« 2 



(24) 
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Fig. 3. Influence of the tail on the cooling coefficient. The symbols show the same data as in Fig. 
[2] (but in linear scale). The full line shows 2(70 — 7) as given by Eqs. (1 1 8 I t and 1261 . 



Performing the integration in Eq. (|23p we obtain 

K -l 42 t(0) 1 ^ 2AB 8 ^ T(k ■ 4 bc ) ( 25) 

where T(a;, y) is the incomplete gamma-function. Inserting into Eq. (|16p yields the 
cooling rate, 

= (1-e 2 ) (1 + ^a 2 ) Ab 6 + 27rV2gr(6, be,) 
7 6(1- I tt2 ) Ab e + 8^V2B6 2 r(4, 6c,) ' 



Naturally, for B = 0, corresponding to absence of the tail (see Eq. (fT0|) ). Eq. (|26|) 
reduces to the previous relation Eq. (fl8|) . In Fig. [3] the difference (70 — 7) which 
quantifies the contribution to the cooling coefficient from the tail is compared with 
the numerical results. 

The scaling law, 70 — 7dsmc exp(— 6e), shown in Fig.[2]as a numerical result is, 
however, difficult to confirm analytically. The approximate analytical theory shown 
in Fig. [3] agrees with the numerical results only qualitatively. 



4. Slow relaxation of the velocity distribution 

In the above analysis we assumed that the velocity distribution function /(c) had 
already relaxed to its steady state scaling shape. In our numerical experiments 
we observed, however, that the relaxation to the stationary form occurs extremely 
slowly, as compared to the relaxation of molecular gases to equilibrium distribu- 
tion. To study this retarded relaxation quantitatively, we use the coefficient 7dsmc 
described above, and define the temperature Tgt(r) oc exp(— 27dsmc t). By defini- 
tion, for t>1 one obtains Tdsmc ~ ^fit since 7dsmc was determined as the best 
exponential fit to Tdsmc(t) for t ^> 1. Hence, the quantity 1 — Tdsmc(t) /7flt(r) 
characterizes the relaxation of the distribution function to its stationary form. We 
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observe that the relaxation time depends sensitively on the coefficient of restitution 
e. 

Figures [4j7] illustrate the relaxation kinetics for different values of the coefficient 
of restitution e. In Figs. [4] and [5] the initial relaxation is shown. We observe a initial 
quick relaxation with characteristic time of 3-5 collisions per particles, similar as 
for molecular gases. 

For values e < 0.7 the simulated temperature approaches Haff's law from below 
(Fig. |4|), whereas for e > 0.7 it approaches Haff's law from above (note the different 
labeling of the vertical axis in these figures). For e = 0.7 the initial relaxation rate 
vanishes. This may be explained by the fact that at e « 0.7 the value of the second 



100 



10 









\. £=0.1 
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e=0.3 
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\ £=0.A 

IN, l V 
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Fig. 4. Short-time relaxation of the temperature decay to Haff's law, Eq. J4) for a system of 
N = 10 s particles for different coefficients of restitution, e = 0.1, e = 0.3, e = 0.4, and e = 0.6. 
The initial relaxation occurs on the time scale 1 — TbsMc/Tflt °^ exp(— 8/15r). 




Fig. 5. Same as Fig. [4] but for e = 0.7, e = 0.8, and e = 0.9. Note that the vertical axis has 
opposite sign as the axis in Fig. [4] 
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Sonine coefficient changes its sign, see Eq. That is, for e < 0.7 the stationary 
velocity distribution / is bent towards lower velocities as compared with the Maxwell 
distribution while for e > 0.7 the distribution is bent towards higher velocities. For 
e = 0.7 the second Sonine coefficient is very small, 02 (0.7) w 0, therefore, for c ~ 1 
the distribution function is very close to the Maxwell distribution and, thus, there 
is no initial relaxation. These arguments prove that the initial relaxation shown 
in Figs. |4] and [5] corresponds to the relaxation of the main part of the velocity 
distribution, c ~ 1, whose deviation from the Maxwell distribution is described by 
the second Sonine coefficient a-i. 

For small coefficient of restitution, e < 0.6, the initial slope of the relaxation 
curves is almost independent of e. Its value is surprisingly close to the slope —8/15 
of the relaxatio n cu rve of the second Sonine coefficient in the limit of almost elastic 
particles, e < Also surprisingly, for larger values of the restitution coefficients 
e > 0.7, which correspond to negative values of the second Sonine coefficient, the 
slope becomes smaller and depends noticeably on e. For almost elastic particles, 
e = 0.9 we find the initial slope 1/2. Presently we do not have an explanation for 
this behavior of the initial relaxation. 

In Figs. [6] and [7] we present the complete relaxation to the steady state. The 
complete relaxation requires much longer time as compared to molecular gases. This 
is related to the relatively slow formation of the exponential tail. The relaxation time 
depends on the coefficient of restitution as shown in Fig. [6j The larger the coefficient 
of restitution the longer it takes to form the complete velocity distribution including 
the exponential tail. 

To explain this effect we use the e quatio n for the relaxation of the velocity 
distribution function to its scaling formj ' 

f (3 + c^j f(c, t) + Jo^Rc, t) = I (/, /) (27) 




Fig. 6. Long-time relaxation of the temperature decay rate for e = 0.1, e = 0.3, and e = 0.4. 
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where /(/,/) is the reduced collision integral, Eq. ([9]), and Jo is given in Eq. (|17p . 
As usual we neglect the incoming term for c 3> lSH leading to the approximation 



I 



(/,/) «-7rc/(c) 



for 



c> 1 



With the Ansatz /(c, r) = _Bexp [— w(t)c] we recast Eq. (|2~T1) into 

dw fj,2 



dr 3Jr 



■ w 



TV 



for 



c > 1 



with the solution 



w(t) = b + (1 — b) exp 



ro(e) 



(28) 



(29) 



(30) 



where b = 3tt//12 coincides with Eq. ([7]) and Tq 1 (e) = ^2/3 Jo- Neglecting 02, which 
quantifies (small) deviations of the main part of the distribution with respect to the 
Maxwellian, and the contribution from the tail, Eq. (pM)) yields Jo = 2^/2n, which 
together with Eq. ([8|) leads to the relaxation time 

-2 



1 



6 



(31) 



For example, for s = 0.4 we obtain the theoretical relaxation time, To ~ 7.1. As 
shown in Fig.[6]for e = 0.4 in the simulation the quantity (1 — TbsMc/Jfit) decreases 
by the factor of 10 in the time span At = 25, ranging from r = 10 to r = 35. This 
gives the estimate for the relaxation time To = 25/log(10) ~ 10.8, in agreement 
with the above theoretical value t = 7.1. 

From the above theory we expect that the relaxation time increases with e. 
While this tendency is confirmed for small coefficients of restitution from e = 0.1 to 
e = 0.4, Fig. [51 the relaxation time seems to saturate for larger e > 0.6, Fig. [7] This 
is, presumably, a finite size effect: For large values of e the number of particles in 
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Fig. 7. Same as Fig.[6]but for e = 0.6, e = 0.7, e = 0.8, and e = 0.9. 
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the tail, which is determined by the threshold velocity c„, increasing with e, is not 
large enough to develop a significant tail. Therefore for such systems the apparent 
relaxation occurs much faster then it would be in a sufficiently large system. The 
relaxation time is mainly determined by the number of particles in the tail, rather 
then by the coefficients of restitution e. Nevertheless, even for these systems, which 
are not sufficiently large for a quantitative study of the tail relaxation, the latter 
occurs anomalously slow. 

5. Conclusion 

We studied analytically and numerically the impact of the high-energy tails of the 
velocity distribution function in granular gases in the homogeneous cooling state on 
the cooling coefficient and on the relaxation time towards the steady state velocity 
distribution. 

In our analytical theory we used an universal Ansatz for the velocity distribution 
function for the entire range of velocities. This Ansatz comprises the main part of the 
distribution function, where the velocities are of the order of the thermal velocity, 
v <~ vt as well as the tail, where v » vt- 

We derived the coefficients of the proposed Ansatz, which allows to estimate 
the cooling coefficient and to characterize the impact of the high-energy tail on the 
cooling. Our analytical results are in qualitative agreement with numerical data 
obtained by Direct Simulation Monte Carlo of 10 8 particles. 

In the simulations we found anomalously slow relaxation of the velocity distri- 
bution function to its stationary form as compared to the relaxation of molecular 
gases, where the relaxation occurs during few collision times. Instead, for granular 
gases, we observed much longer relaxation, in the order of 20-30 collisions per par- 
ticle. To explain this slow relaxation we developed a theory, which predicts that the 
relaxation time increases with decreasing inelasticity, in agreement with the numer- 
ical observations. For a large coefficient of restitution the relaxation time saturates 
with increasing e, which may be attributed to finite size effects, that is, the system 
of 10 s particles is not large enough for the quantitative numerical analysis of the 
tail relaxation. 
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